Climate: Land Use

Introduction

This script contains the methods used to wrangle, analyse and present land use data in the Northern Three regions. For a guide on downloading land use data refer to the README document for the Spatial Analysis GitHub repo. Note that the data for this script requires an additional step of pre-processing using QGIS.

Land use data is predominantly used within the climate section of the technical report to “set the scene” for each basin in each region. Land use is currently not scored, and is only a contextual dataset. The main objectives of this script are to:

• Create tabular summaries of the types of land use present in each basin.
• Create a plot of different land use types as a proportion of the basin.
• Create maps of land use types in each basin.

Script Set Up

This script requires multiple core spatial packages, each of which is loaded below. Key variables and paths are also created.

#use pacman function to load and install (if required) all other packages
pacman::p_load(dplyr, tidyr, readr, purrr, stars, glue, here, sf, tmap, RColorBrewer, ggplot2, stringr, lubridate)

#install/load the custom RcTools package
#pak::pak("add-am/RcTools")
library(RcTools)

#define the script crs and target year
script_crs <- "EPSG:7855"
script_fyear <- 2025
script_fac <- 5

sf_use_s2(FALSE)

#build the data path and output path variables
data_path <- here("data/land_use/")
output_path <- here(glue("outputs/land_use/{script_fyear}/"))

#and create some additional sub folders
dir.create(glue("{data_path}/raw/"), recursive = TRUE)
dir.create(glue("{data_path}/processed/"))
dir.create(glue("{output_path}/plots/"), recursive = TRUE)
dir.create(glue("{output_path}/maps/"))

Load Data

Now the script is set up we need to load in all of the required datasets. This will be broken into two segments: - Spatial data specific to the N3 region - such as the region, basin, and sub basin boundaries. - Land use data

Spatial Data

The northern three boundaries are built using a custom function or reloaded from file.

#if the data is on file, load it, otherwise built it and save it for next time
if (file.exists(here("data/n3_region.gpkg"))){
  n3_region <- st_read(here("data/n3_region.gpkg"))
} else {
  n3_region <- build_n3_region()
  st_write(n3_region, here("data/n3_region.gpkg"))
}

The data is then filtered for ease of use later.

#select only the land components of the n3 region dataset
n3_sub_basins <- n3_region |> 
  filter(Environment != "Marine") |> 
  group_by(Region, BasinOrZone, SubBasinOrSubZone) |> 
  summarise(geometry = st_union(geom)) |> 
  ungroup() |> st_cast()

Land Use

For the N3 region we require multiple large datasets (e.g. across the entirety of Queensland) that can take a while to process. The first time these code chunks are run they will read in each of the land use datasets, perform the editing steps, and then save the output. The next time the code chunks are run they will check if the edited version exists, and open that instead. This dramatically reduces processing time on repeat runs of the script.

Note

If for some reason you find that the edited output is not what you want, or that you need to re-run the original code chunks you will need to edit the code chunk to not allow the edited version to be opened. Or delete the edited version from the directory.

For example, if new datasets are published to QSpatial you will need to rerun the original merge and crop step.

All land use data is downloaded from QSpatial:

• Search for:
  • “Land use mapping - 1999 to 2016 - Burdekin NRM”.
  • “Land use mapping - 1999 to 2015 - Wet Tropics NRM region”.
  • “Land use mapping - 1999 to 2016 - Mackay Whitsunday NRM”.
  • “Land use mapping - 2021 - Great Barrier Reef NRM regions”.
• Download as is, there are no formatting options.
• Once downloaded, for the Burdekin, Wet Tropics, and GBR datasets, extract the entire .gdb folder and rename each .gdb folder to the following:
  • “dt_land_use.gdb”
  • “wt_land_use.gdb”
  • “qld_land_use.gdb”
• For the Mackay Whitsunday NRM*, extract the .gpkg file and rename the file to “mwi_land_use.gpkg”.

A custom function has then been written to handle the extraction and refining of each of these datasets into a single object.

#build path to the raw data
raw_data_path <- glue("{data_path}/raw/")

#build path to where cropped datasets should be saved
cropped_data_path <- glue("{data_path}/processed/")

#run custom function
n3_land <- extract_land_use(raw_data_path, cropped_data_path, n3_sub_basins)

#names need to be made lowercase for simplicity later
n3_land <- n3_land |> 
  mutate(across(where(is.character), tolower))

Once the main dataset has been made we can group the dataset into what is currently accepted as the broadest level of classification as these are the classifications used in the Land Use analysis. This is because, generally we are not concerned with the difference between for example; “National Park” and “Natural Feature Protection” (a tertiary level separation). Therefore we will only breakdown and report on land use at a custom mix of the primary and secondary land use type levels. These levels are as follows:

• “Conservation”
• “Dryland Agriculture”
• “Forestry”
• “Grazing”
• “Irrigated Agriculture”
• “Mining”
• “Urban/Intensive”
• “Water”

#aggregate data
n3_land <- n3_land |> 
  select(-Tertiary) |> 
  mutate(across(where(is.character), tolower)) |> 
  mutate(Landuse = case_when(
    str_detect(Primary, "water") ~ "water",
    str_detect(Primary, "conservation") ~ "conservation",
    str_detect(Primary, "irrigated") ~ "irrigated agriculture",
    str_detect(Primary, "dryland") ~ "dryland agriculture",
    str_detect(Primary, "relatively natural") ~ "grazing",
    str_detect(Primary, "intensive") ~ "urban/intensive")) |> 
  mutate(Landuse = case_when(
    str_detect(Secondary, "grazing") ~ "grazing", 
    str_detect(Secondary, "forest") ~ "forestry",
    str_detect(Secondary, "mining") ~ "mining", 
    TRUE ~ Landuse)) |> 
  select(-c(Primary, Secondary))

Analyse Data

Now the data is sorted, imported, and ready to be reported, we can move onto calculating some statistics. Below we calculate:

• The total area of each region each year (varies very slightly year on year)
  • The area of each landuse type in each region each year
  • The proportion of each landuse type relative to the entire region area each year
  • The area change of each landuse type from year to year (i.e. if 1999 = 2km, and 2009 = 3km, area change = 1km)
  • The percentage change of each landuse type from year to year (i.e. area change = 1km (and 1999 = 2km), percent change = 50%)

The above calculations are then repeated for every basin, and then repeated again for every sub basin.

Note

Currently there is no defined “nice” way that we have elected to present this data, so we are simply saving all the data in one table. However, in the future we may be able to define the format in which we want to present, and exactly what data is needed.

#calculate the area of every single row, then group up and sum areas for each group
n3_land_use_tbl <- n3_land |> 
  mutate(LanduseArea = st_area(geom)) |> 
  st_drop_geometry() |> 
  group_by(Landuse, Region, BasinOrZone, SubBasinOrSubZone, Year) |> 
  summarise(LanduseArea = sum(LanduseArea)) |> 
  ungroup()

#update the units to a more reasonable metric
n3_land_use_tbl$LanduseAreaKm2 <- units::set_units(n3_land_use_tbl$LanduseArea, km^2)

#drop the original area (m)
n3_land_use_tbl <- n3_land_use_tbl |> dplyr::select(-any_of("LanduseArea"))

#create three different vectors to be used as a grouping and a similar list for removing columns
cols_to_pick <- list(c("Region"), c("Region", "BasinOrZone"), c("Region", "BasinOrZone", "SubBasinOrSubZone"))
cols_to_remove <- list(c("BasinOrZone", "SubBasinOrSubZone"), c("SubBasinOrSubZone"), c(""))

#create an empty df to store the output
n3_land_use_final <- data.frame()

#for each of the grouping vectors
for (i in 1:length(cols_to_pick)){
  
  #calculate everything (i.e. do the math at a region level, a basin level, and a sub basin level)
  temporary_table <- n3_land_use_tbl |> 
    group_by(across(.cols = c(cols_to_pick[[i]], Year))) |> #group at year + region, then next loop at year + basin, then next loop, etc.
    mutate(TotalAreaKm2 = sum(LanduseAreaKm2)) |> #get the total area for the group
    group_by(across(.cols = c(cols_to_pick[[i]], Landuse, Year))) |> #do the same group, including landuse type
    mutate(LanduseKm2 = sum(LanduseAreaKm2), #get the area of the year + landuse + region/basin/sub basin
           LanduseProportionOfTotalPercent = (LanduseKm2/TotalAreaKm2)*100) |> #get the proportion of the landuse vs the total area
    dplyr::select(-any_of(c(cols_to_remove[[i]], "LanduseAreaKm2"))) |> #remove cols not used in the original grouping
    distinct() |> 
    group_by(across(.cols = c(cols_to_pick[[i]], Landuse))) |> #group at land use + region/basin/or sub basin (ignore year)
    mutate(across(LanduseKm2, as.numeric), #force area to be a numerical value for calculations below
           AreaChangeYearToYear = (LanduseKm2 - lag(LanduseKm2)), #calculate the difference between each row (year) as km2
           PercentageChangeYearToYear = (LanduseKm2/lag(LanduseKm2) -1) *100, #calculate the difference as %
           XyearToXyear = paste(Year, lag(Year), sep = " to "), #create column explain what years were compared
           XyearToXyear = case_when(str_detect(XyearToXyear, "NA") ~ "1999 to 1999", #edit to fix when there is no "row above"
                                      T ~ XyearToXyear),
           across(where(is.numeric), \(x) replace_na(x, 0)), #replace na area and % change values with zero
           across(where(is.numeric), \(x) round(x, 1))) |> ungroup()
  
  #join to main
  n3_land_use_final <- bind_rows(n3_land_use_final, temporary_table)

}

#reorder cols
n3_land_use_final <- n3_land_use_final |> 
  relocate(Region, BasinOrZone, SubBasinOrSubZone)

#save the output
write_csv(n3_land_use_final, glue("{output_path}/land_use_full_table.csv"))

Visualise Land Use

We can now visualise land use, using both plots and graphs.

Mapping Land Use

Currently we are presenting land use maps for:

• Each sub basin in each region: DT(12), Burdekin (6), –> 18
• Each basin in each region: DT(2), WT(9), MWI(5), Burdekin (2), –> 18
• Each region overall (DT, WT, MWI, Burd), –> 4
• and, Across every year: (which is currently 1999, 2009, 2015/2016 and 2021). –> 4

for a total of (18+18+4)*4 = 160 maps.

First we will prep the data and make sure all the correct colours are assigned.

#create a named colour palette
my_palette <- c("conservation" = "#277402", 
                "dryland agriculture" = "#948A54", 
                "forestry" = "#4EE703",
                "grazing" = "#FFFFBF", 
                "irrigated agriculture" = "#D0FF73", 
                "mining" = "#632523",
                "urban/intensive" = "#FF5704", 
                "water" = "#4F81BD")

#update the dataset by calling the palette colours using the landuse column in the df
n3_land <- n3_land |> 
  mutate(Palette = my_palette[Landuse])

Then we will create a list contains name vectors of targets. Essentially providing us with the name of the column to filter, and the name of the location to look for.

#get a table of unique locations that need to be mapped
unique_locations <- n3_land |> 
  st_drop_geometry() |> 
  select(Region, BasinOrZone) |> 
  unique()

#build a list of basins within each Region
basins_per_region <- as.list(split(unique_locations$BasinOrZone, unique_locations$Region))

#join on the the list of just basins
all_locations <- c(basins_per_region, unique_locations$BasinOrZone) 

#name each item in the list
names(all_locations) <- c(names(basins_per_region), unique_locations$BasinOrZone)

#and edit the one duplicate in the naming
names(all_locations)[1] <- "burdekin_region"

Using this list of targets we can loop over the dataset to create each map.

#map over the unique locations and create region maps
purrr::walk2(all_locations, names(all_locations), \(x, y) {

  #filter the main dataset
  target_data <- n3_land |> 
    filter(BasinOrZone %in% x)

  #create an outline of the same location
  target_outline <- target_data |> 
    filter(Year == min(Year)) |> 
    st_union() |> 
    st_as_sf() |> 
    st_collection_extract("POLYGON")

  #plus an outline of the broader region
  target_region <- n3_land |> 
    filter(
      Region == unique(target_data$Region),
      Year == min(Year)) |> 
    st_union() |> 
    st_as_sf() |> 
    st_collection_extract("POLYGON")

  #create the water layer for the location
  water_layer <- maps_water_layer(x, WaterLakes = TRUE, StreamOrder = 3)

  #create the inset objects for the location
  inset_objects <- maps_inset_layer(target_outline, target_region, 1.1)

  #create a string defining the save folder
  save_folder <- glue("{output_path}/maps/{unique(target_data$Region)}/")

  #create the folder if it doesn't already exist
  if(!dir.exists(save_folder)){dir.create(save_folder)}

  #for each year present in the filtered data, create a map
  for (i in unique(target_data$Year)){

    #select a year
    single_year <- filter(target_data, Year == i)

    #create the map
    lu_map <- tm_shape(qld) +
      tm_polygons(fill = "grey80", col = "black") +
      tm_shape(target_region) +
      tm_polygons(fill = "grey90", col = "black") +
      tm_shape(single_year, is.main = TRUE) +
      tm_fill(fill = "Palette") +
      tm_add_legend(
        type = "polygons", 
        fill = unique(single_year$Palette), 
        labels = unique(single_year$Landuse),
        title = "Land Use") +
      water_layer +
      tm_shape(target_outline) +
      tm_borders(col = "black") +
      tm_layout(
        legend.frame = TRUE, 
        legend.bg.color = "White", 
        legend.text.size = 0.7, 
        legend.position = c("left", "bottom"),
        asp = 1.1)

    #save the map
    tmap_save(
      lu_map, 
      filename = glue("{save_folder}/{y}_{i}_land_use.png"),
      insets_tm = inset_objects[1], 
      insets_vp = inset_objects[2])
  }
})

As an example of what one of these maps looks like, see below:

map

Mapping Land Use Change

Generally, it is equally important to understand not only what the current or historical landuse is, but also how the landuse has changed over time. A custom function has been written that conducts this analysis

#compare two years of data
change_in_use <- sf_change_over_time(n3_land, "Landuse", 1999, 2021)

#intsect back over original to retrieve important metadata information
change_in_use <- st_intersection(change_in_use, n3_sub_basins) |> 
  mutate(across(where(is.character), tolower))

In the case of landuse, we end up with 64 unique “change possibilities”, (which are the cross cases of the 8 land uses times the 8 land uses (e.g. Grazing to Conservation)). Although useful in tabular form, trying to create a map with 64 unique classifications is not the best idea, so below we create some maps focusing on the changes that most interest us, specifically for:

• Urban/Intensive
  • divided more specifically into:
    • Urban to urban (colour = W) (No Change to Urban)
    • Urban to anything not urban (colour = X) (Loss of Urban)
    • Anything not urban to urban (colour = Y) (Gain of Urban)
    • Anything that is not urban in either dataset (colour = Z) (Not Urban)
• Conservation
  • divided more specifically into:
    • Conservation to conservation (colour = W) (No Change to conservation)
    • Conservation to anything not conservation (colour = X) (Loss of conservation)
    • Anything not conservation to conservation (colour = Y) (Gain of conservation)
    • Anything that is not conservation in either dataset (colour = Z) (Not conservation)

#create an urban change focus table and assign colours
urban_use <- change_in_use |> 
      mutate(
        "Change" = case_when(
          str_detect(ValueChange, "urban/intensive to urban/intensive") ~ "No Change",
          str_detect(ValueChange, "uUrban/intensive to (?!urban/intensive\\b).*") ~ "Loss",
          str_detect(ValueChange, "(?<!urban/intensive) to urban/intensive") ~ "Gain",
          TRUE ~ "Other Use"),
        "Palette" = case_when(
          Change == "No Change" ~ "#FFFF00",
          Change == "Loss" ~ "#277402",
          Change == "Gain" ~ "#FF0000",
          Change == "Other Use" ~ "#A6A6A6"))

#create a conservation focused table and assign colours
conservation_use <- change_in_use |> 
      mutate(
        "Change" = case_when(
          str_detect(ValueChange, "conservation to conservation") ~ "No Change",
          str_detect(ValueChange, "conservation to (?!conservation\\b).*") ~ "Loss",
          str_detect(ValueChange, "(?<!conservation) to conservation") ~ "Gain",
          TRUE ~ "Other Use"),
        "Palette" = case_when(
          Change == "No Change" ~ "#FFFF00",
          Change == "Loss" ~ "#FF0000",
          Change == "Gain" ~ "#277402",
          Change == "Other Use" ~ "#A6A6A6"))

#combine the two datasets into a list
focused_landuse <- list(urban_use, conservation_use)

#and give them names
names(focused_landuse) <- c("urban", "conservation")

#map over the two dataset types, the (nested) map over over the unique locations and create maps
purrr::walk2(focused_landuse, names(focused_landuse), \(df, df_name){
  purrr::walk2(all_locations, names(all_locations), \(location, location_name) {

    #filter the focused landuse dataset
    target_data <- df |> 
      filter(BasinOrZone %in% location)

    #create an outline of the same location
    target_outline <- target_data |> 
      st_union() |> 
      st_as_sf() |> 
      st_collection_extract("POLYGON")

    #plus an outline of the broader region
    target_region <- n3_land |> 
      filter(
        Region == unique(target_data$Region),
        Year == min(Year)) |> 
      st_union() |> 
      st_as_sf() |> 
      st_collection_extract("POLYGON")

    #create the water layer for the location
    water_layer <- maps_water_layer(location, WaterLakes = TRUE, StreamOrder = 3)

    #create the inset objects for the location
    inset_objects <- maps_inset_layer(target_outline, target_region, 1.1)

    #create a string defining the save folder
    save_folder <- glue("{output_path}/maps/{unique(target_data$Region)}/")

    #create the folder if it doesn't already exist
    if(!dir.exists(save_folder)){dir.create(save_folder)}

    #create the map (dont forget to add the custom water layer)
    lu_change_map <- tm_shape(qld) +
      tm_polygons(fill = "grey80", col = "black") +
      tm_shape(target_region) +
      tm_polygons(fill = "grey90", col = "black") +
      tm_shape(target_data, is.main = TRUE) +
      tm_fill(fill = "Palette") +
      tm_add_legend(
        type = "polygons", 
        fill = unique(target_data$Palette), 
        labels = unique(target_data$Change),
        title = "Landuse Change") +
      water_layer +
      tm_shape(target_outline) +
      tm_borders(col = "black") +
      tm_layout(
        legend.frame = TRUE, 
        legend.bg.color = "White", 
        legend.text.size = 0.7, 
        legend.position = c("left", "bottom"),
        asp = 1.1)
        
    #save the map as a png
    tmap_save(
      lu_change_map, 
      filename = glue("{save_folder}/{location_name}_{df_name}_use_change.png"),
      insets_tm = inset_objects[1], 
      insets_vp = inset_objects[2])

  })
})

Plot Land Use

The other visualisation that is useful to gain perspective on land use in each basin is to plot the most up to date land use in a stacked bar chart.

#use the dataset built during the initial analysis phase
n3_land_use_final <- n3_land_use_final |> 
  mutate(
    LanduseProportionOfTotalPercent = as.numeric(LanduseProportionOfTotalPercent), 
    Year = as.numeric(Year))

#assign names to my palette so ggplot can colour everything correctly
names(my_palette) <- unique(n3_land_use_final$Landuse)

for (i in unique_locations$BasinOrZone){
  
  #select a specific basin and the latest year
  basin_plot <- n3_land_use_final |> 
    filter(BasinOrZone == i, Year == max(n3_land_use_final$Year)) |> 
    arrange(desc(LanduseProportionOfTotalPercent)) |> 
    mutate(Landuse = factor(Landuse, levels = unique(Landuse)))
  
  #create a plot
  plot <- ggplot(basin_plot, aes(fill = Landuse, x = Year, y = LanduseProportionOfTotalPercent)) +
    geom_bar(position = position_fill(reverse = TRUE), stat = "identity") +
    scale_fill_manual(values = my_palette) +
    scale_x_discrete(expand = c(0.025, 0)) +
    scale_y_continuous(expand = c(0, 0), labels = scales::percent_format(accuracy = 1)) +
    labs(x = "", y = "", fill = glue("{i} Basin Land Use")) +
    theme(
      axis.text.x = element_text(colour = "black"),
      axis.text.y = element_blank(),
      axis.line.x = element_line(colour = "black"),
      axis.ticks.length = unit(-0.15, "cm"),
      panel.background = element_blank()) +
    coord_flip()
  
  #save
  ggsave(glue("{output_path}/plots/{i}_land-use.png"), plot, width = 12, height = 4)
}

After running and saving all of these see below for an example.

plot